n_distinct(result$MA.id) # number of MAs
## [1] 23492
nrow(result) # 9 model result per each MA
## [1] 211428
table((result %>% distinct(MA.id, .keep_all = TRUE))$event)
##
## common infrequent rare very rare
## 22074 1154 101 163
proportions(table((result %>% distinct(MA.id, .keep_all = TRUE))$event)) * 100
##
## common infrequent rare very rare
## 93.9639026 4.9123106 0.4299336 0.6938532
There are three types of model that have \(\tau\)
ivre,
htevaguelogit and htelogit_tau50 (log
ratios)ablogit and
ablogit_tau50 (log odds)htebeta
(prob)cor(
result_bayesian[result_bayesian$modelname == "htelogit_tau50",]$tau,
result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$taua
)
## [1] 0.2477361
cor(
result_bayesian[result_bayesian$modelname == "htelogit_tau50",]$tau,
result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$tauc
)
## [1] 0.3046824
cor(
result_bayesian[result_bayesian$modelname == "htebeta",]$taua,
result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$taua
)
## [1] -0.08585503
cor(
result_bayesian[result_bayesian$modelname == "htebeta",]$tauc,
result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$tauc
)
## [1] -0.1115232
| modelname | Sig_sum | percentage_total |
|---|---|---|
| ivfe | 10403 | 44.283160 |
| ivre | 9132 | 38.872808 |
| ctevaguelogit | 11553 | 49.178444 |
| ctebeta | 11402 | 48.535672 |
| htevaguelogit | 7308 | 31.108462 |
| htelogit_tau50 | 6047 | 25.740678 |
| htebeta | 4363 | 18.572280 |
| ablogit | 8003 | 34.066916 |
| ablogit_tau50 | 597 | 2.541291 |
table(result_sigsum$Sig_sum)
##
## 0 1 2 3 4 5 6 7 8 9
## 10947 1075 706 1407 1587 922 915 2312 3121 500
round(proportions(table(result_sigsum$Sig_sum))*100,3)
##
## 0 1 2 3 4 5 6 7 8 9
## 46.599 4.576 3.005 5.989 6.755 3.925 3.895 9.842 13.285 2.128
table(result_sigsum$event, result_sigsum$Sig_sum)
##
## 0 1 2 3 4 5 6 7 8 9
## common 9948 977 632 1340 1524 902 883 2261 3109 498
## infrequent 799 74 58 59 61 18 29 42 12 2
## rare 65 8 7 5 2 2 3 9 0 0
## very rare 135 16 9 3 0 0 0 0 0 0
round(proportions(table(result_sigsum$event, result_sigsum$Sig_sum), margin = 1)*100,2)
##
## 0 1 2 3 4 5 6 7 8 9
## common 45.07 4.43 2.86 6.07 6.90 4.09 4.00 10.24 14.08 2.26
## infrequent 69.24 6.41 5.03 5.11 5.29 1.56 2.51 3.64 1.04 0.17
## rare 64.36 7.92 6.93 4.95 1.98 1.98 2.97 8.91 0.00 0.00
## very rare 82.82 9.82 5.52 1.84 0.00 0.00 0.00 0.00 0.00 0.00
For best model comparison among seven Bayesian models (exclude ivfe and ivre)
result_htevaguelogit_tau = result %>% filter(modelname %in% c("ctevaguelogit", "htevaguelogit")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>% filter(modelname == "htevaguelogit") %>% ungroup()
mean(result_htevaguelogit_tau$diff_waic < 0) # htevaguelogit has better waic than ctevaguelogit
## [1] 0.9258471
result_htelogit_tau50_tau = result %>% filter(modelname %in% c("ctevaguelogit", "htelogit_tau50")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>% filter(modelname == "htelogit_tau50") %>% ungroup()
mean(result_htelogit_tau50_tau$diff_waic < 0) # htelogit_tau50 has better waic than ctevaguelogit
## [1] 0.919164
result_htelogit = result %>% filter(modelname %in% c("htevaguelogit", "htelogit_tau50")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>% filter(modelname == "htelogit_tau50") %>% ungroup()
mean(result_htelogit$diff_waic < 0) # htelogit_tau50 has better waic than htevaguelogit
## [1] 0.5440576
result_ablogit_tau = result_narm %>% filter(modelname %in% c("ctevaguelogit", "ablogit")) %>% group_by(MA.id) %>%
mutate(diff_waic = diff(waic)) %>% filter(modelname == "ablogit") %>% ungroup()
mean(result_ablogit_tau$diff_waic < 0) # ablogit has better waic than ctevaguelogit
## [1] 0.5141309
result_ablogit_tau50_tau = result_narm %>%
filter(modelname %in% c("ctevaguelogit", "ablogit_tau50")) %>% group_by(MA.id) %>%
mutate(diff_waic = diff(waic)) %>% filter(modelname == "ablogit_tau50") %>% ungroup()
mean(result_ablogit_tau50_tau$diff_waic < 0) # ablogit_tau50 has better waic than ctevaguelogit
## [1] 0.6156923
result_ablogit_tau = result_narm %>% filter(modelname %in% c("ctebeta", "ablogit")) %>% group_by(MA.id) %>%
mutate(diff_waic = diff(waic)) %>% filter(modelname == "ablogit") %>% ungroup()
mean(result_ablogit_tau$diff_waic < 0) # ablogit has better waic than ctebeta
## [1] 0.7459887
result_ablogit_tau50_tau = result_narm %>%
filter(modelname %in% c("ctebeta", "ablogit_tau50")) %>% group_by(MA.id) %>%
mutate(diff_waic = diff(waic)) %>% filter(modelname == "ablogit_tau50") %>% ungroup()
mean(result_ablogit_tau50_tau$diff_waic < 0) # ablogit_tau50 has better waic than ctebeta
## [1] 0.7681851
result_ablogit = result_narm %>%
filter(modelname %in% c("ablogit", "ablogit_tau50")) %>% group_by(MA.id) %>%
mutate(diff_waic = diff(waic)) %>% filter(modelname == "ablogit_tau50") %>% ungroup()
mean(result_ablogit$diff_waic < 0) # ablogit_tau50 has better waic than ablogit
## [1] 0.6549662
## ctevaguelogit ctebeta htevaguelogit ablogit htebeta
## ctevaguelogit 0.000 0.741 0.074 0.486 0.280
## ctebeta 0.259 0.000 0.062 0.254 0.040
## htevaguelogit 0.926 0.938 0.000 0.920 0.754
## ablogit 0.514 0.746 0.080 0.000 0.265
## htebeta 0.720 0.960 0.246 0.735 0.000
## htelogit_tau50 0.919 0.945 0.541 0.925 0.782
## ablogit_tau50 0.616 0.768 0.199 0.655 0.443
## htelogit_tau50 ablogit_tau50
## ctevaguelogit 0.081 0.384
## ctebeta 0.055 0.232
## htevaguelogit 0.459 0.801
## ablogit 0.075 0.345
## htebeta 0.218 0.557
## htelogit_tau50 0.000 0.828
## ablogit_tau50 0.172 0.000
## ctevaguelogit ctebeta htevaguelogit ablogit htebeta
## ctevaguelogit 0.000 0.000 0.000 0.000 0.000
## ctebeta 0.259 0.000 0.000 0.000 0.000
## htevaguelogit 0.926 0.938 0.000 0.000 0.000
## ablogit 0.514 0.746 0.080 0.000 0.000
## htebeta 0.720 0.960 0.246 0.735 0.000
## htelogit_tau50 0.919 0.945 0.541 0.925 0.782
## ablogit_tau50 0.616 0.768 0.199 0.655 0.443
## htelogit_tau50 ablogit_tau50
## ctevaguelogit 0.000 0
## ctebeta 0.000 0
## htevaguelogit 0.000 0
## ablogit 0.000 0
## htebeta 0.000 0
## htelogit_tau50 0.000 0
## ablogit_tau50 0.172 0